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INTRODUCTION 


PROBLEM  DESCRIPTION 

Consider  the  motion  of  a body  (or  point  mass)  in  an  inverse- square 
gravitational  field  centered  at  the  origin  of  some  arbitrary  inertial  ref- 
erence frame.  At  any  point  R the  body  will  experience  an  acceleration  of 

2 “ 3 2 

-K  R R,  where  K is  the  gravitational  constant  for  the  given  field. 

-f 

Now  suppose  two  nonparallel  position  vectors  R,  R>p  and  a time  tf 
are  given.  The  correlated  velocity  (denoted  for  these  conditions  is 

-4  4 

defined  as  the  velocity  at  R which  will  cause  a body  to  reach  R ip  in  time 

4 -4  -4 

tf  if  no  other  force  acts  on  it*  Clearly  Vr  is  determined  by  R,  tf/ 

2 4-4  -4  -4  « -*  4-4,. 

and  K , so  we  may  write  Vc  = Vc(tf,  R,  Rip,  K ).  The  symbols  3vc/3r,  3vc/3tf, 
etc.,  will  be  used  to  denote  the  partials  obtained  by  holding  the  other  pa- 

-4  -4 

rameters  in  this  set  constant.  In  particular,  dV^/dR  is  commonly  called 
the  Q-matrix,  and  denoted  simply  by  Q. 

It  is  well  known  that  the  path  of  a body  in  any  inverse- square  force  field 
will  lie  along  a conic  section  with  the  center  of  force  as  a focus.  We  will 
consider  only  the  case  in  which  the  path  is  along  an  ellipse  and  subtends  an 

-4 

angle  less  than  tt.  vc  will  be  considered  undefined  if  no  value  corresponding 

-4 

to  such  a path  exists.  This  restriction  assures  that  Vc  is  well-defined 
when  it  exists. 

This  paper  describes  a new  method  for  computing  the  Q-matrix  which  is 
completely  analytic,  yet  is  also  simple  and  efficient. 

BASIC  CONCEPTS  AND  TERMINOLOGY 

-4  -4 

Clearly  R and  R^,  determine  the  plane  of  the  trajectory.  It  is  con- 
venient to  define  an  orthogonal  coordinate  system  (which  will  be  called  the 
local  coordinate  frame)  with  the  axes 


U j = unit  (R), 


= unit  [ (R  x Rp)  x R] , and 


U3  = X V 


L 


We  then  define  V „ = U,  • V and  V „ = U. 

CR  1 C CO  2 


vc  = °) . 


Similarly,  VR  = and  V^=  U.,  • V,  where  V is  the  actual  velocity 

at  R.  The  range  angle,  0R,  is  the  angle  between  R and  R^.  Finally,  the 
difference  vc  ~ V will  be  called  the  velocity  to  be  gained  and  denoted  by  Vg. 

Since  VCR  and  Vcg  are  scalars,  it  is  often  desirable  to  express  them 

as  functions  of  scalars.  One  such  set  is  (tf,  R,  R^,  0R,  K2);  unless 

otherwise  indicated  V and  V will  be  treated  as  functions  of  these  param- 

cr  cy 

eters.  Thus  9v  „/9tf,  9v  /90  , etc.,  will  denote  the  partials  obtained  by 
CU  i CR  K 

holding  the  other  parameters  in  this  set  constant. 


The  usual  point  of  view  is  that  R is  the  current  position  and  Rt  the  pos- 
ition at  which  the  "target"  will  be  located  after  an  elapsed  time  t^  (often 
called  the  "time-to-go") . Thus  is  the  velocity  required  to  reach  the  tar- 
get in  the  specified  time.  With  this  interpretation  it  makes  sense 
to  consider  R,  V , VCR#  and  V^g  as  functions  of  time  (t) . In  this  case  it 
is  understood  that  t = 0 at  R and  t^  decreases  with  time;  i.e.,  dt^/dt  = -1. 
(The  idea  is  that  we  want  to  reach  the  target  at  a specified  time.)  Of 
course  and  V^R  are  always  defined  relative  to  the  current  position;  that 

is,  V (t)  = V (t)  • unit[R(t)],  and  similarly  for  V ».  Unless  otherwise 
CR  C ->■-*.  ->. 

indicated,  it  is  assumed  that  V = V . Note  that,  if  V(0)  = V (0)  and  no 
force  but  gravity  is  acting  on  the  body,  then  V(t)  = V^ft)  for  0 < t < t^. 

APPLICATIONS 

The  algorithm  developed  in  this  paper  is  useful  in  the  analysis  of 
guidance  systems  using  V steering.  In  particular,  it  can  be  used  to  compute 
the  time  derivative  of  Vg  by  means  of  equation  Bl,  which  is  derived  in  Ap- 
pendix B. 

The  equations  for  9vcg/9tf  and  ^Vcg/90R  are  also  of  interest  in  their 
own  right. 


The  algorithm  for  computing  Q and  the  equation  for  9V(_,q/90r  are  used 
in  several  trajectory  simulation  programs  developed  by  NSWC. 


r- 


mm 


DEVELOPMENT  OF  Q- MATRIX  ALGORITM 


PRELIMINARY  RESULTS 

Vg  steering  and  the  Q-matrix  have  both  been  studied  extensively,  and 
some  earlier  results  will  be  useful  here. 

1.  An  efficient  algorithm  for  computing  t^  for  a given  value  of 
VC0  waS  devel°Ped  by  T.  Alexander.  Differentiation  of  the  equations  involved 
also  yields  an  analytic  expression  for  Sv^/St^.  A derivation  of  these  equa- 
tions, together  with  a suggested  iteration  method  for  implementing  them,  is 
given  in  Appendix  A. 

One  of  the  equations  developed  in  this  derivation,  e.g.  (A5) , can  be 

rearranged  to  give  a useful  formula  for  V : 

CK 


2.  The  Q-matrix  is  symmetric.* 


DERIVATION 


The  algorithm  will  be  developed  in  two  stages.  First,  a representation 
for  Q in  the  local  coordinate  frame  in  terms  of  V^g,  VCR'  and  certain  of 


their  derivatives  will  be  obtained, 
vativer  will  be  derived. 

To  obtain  the  representation  of 

• "*0 
which  Q is  to  be  evaluated.  Let  U , 


Second,  formulas  for  the  required  deri- 
Q,  fix  R and  t^,  and  a position  RQ  at 

— >■  A A 

U“ , and  U”  be  the  axes  of  the  local 


* C.  J.  Cohen  and  R.  H.  Lyddane,  The  Q-Guidance  Matrix  and  Its  Symmetry  (U) , 
Naval  Proving  Ground,  Dahlgren,  Va.,  TR  K-6/59,  June  1959.  UNCLASSIFIED 
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coordinate  system <£  defined  at  RQ , and  let  0po  be  the  corresponding  range 

-V 

angle.  For  an  arbitrary  position  R,  let 
-*•  -*-g 

RR  = R . U°, 

R0  = R • U®,  and 

-+A 

Rp  * R • U3  • 

■>  -V 

For  any  vector  X,  the  representation  of  X in will  be  denoted  by  X , and  the 

y 

corresponding  representation  for  Q by  Q . 

-V  -V  -f 

N.B.:  While  the  axes  Uj»  , and  U are  always  defined  relative  to  R,  the 


symbol denotes  the  fixed  frame  U°,  U°,  U 
By  definition 

8 - 557'^-  aH*  j 

\ R 0 P /R  = R„ 


Also,  it  follows  immediately  from  the  definitions  that  = V^R  Uj  + V 


Applying  the  chain  rule  to  this  equation  yields 


CJ> 

L<+ 

0 

11 

9v 

CR  - 

8V 

+ V + 

9r. 

9r.  ’ 

CR  9r.  9r. 

U2  + VC0 


where  R^  is  any  component  of  R (e.g.,  RR,  R0,  R^) 
It  is  easily  verified  that 


0.  = R- 1 R 


U2  = (Rt  - R~ 1 cos  0R  R)/sin  0 , where  ^ = unit(RTK  ) 

Since  R = (cos  0 , sin  0 , 0)T,  in  X.  these  equations  take  the  form 

T Rq  Rq  . 


R"1  K 


R-1  R. 


(cos  0R  - R-1  RR  cos  0R)csc  8f 


of  =(  R_1  Re  } of  = (sin  9r0  - r_1  r0  cos  eR)csc  eR 


-R-1  R^  cot  0f 


The  same  operation  can  now  be  performed  on  Equation  5 to  give  (at  R0); 


Next  recall  that  Vcg  and  VCR  are  functions  of  tj,  R,  RT,  0R,  and  K2 . Since 
tf,  Rt,  and  K2  are  constant  here.  Equation  6 yields  immediately 

9vci  avci  3vci  _i  3vci  svci 

3r^~  = ~3r~  ; 3r^~  = "R  36^"  '■  3r^  = 0 ( 

where  Vci  can  be  interpreted  as  either  Vcg  or  VCR. 


Finally,  substituting  Equations  7,  8,  and  9 into  Equation  3 with  Ri 
interpreted  as  Rr,  Rq,  and  Rp  and  substituting  the  results  into  Equation  2 
yields 


This  is  the  representation  of  Q mentioned  earlier.  We  now  turn  to  the 

problem  of  finding  formulas  for  the  partials  of  V and  V „ with  respect  to 

CR  C0 

R and  0 . 

R 

First,  Equation  1 gives  V as  a function  of  V , R,  R , 0 , and  K2 , 

CR  Co  T R 

say  V = f (V  R,  R , 0 , K2).  The  derivatives  of  V which  occur  in 
LK  Co  T F.  CR 

Equation  11,  however,  are  defined  with  respect  to  the  set  of  parameters  t^, 
R,  R , 0 , K2 , as  explained  in  the  Introduction.  The  partials  of  f can  be 
obtained  simply  by  differentiating  Equation  1.  This  gives 

Hr  = " VC0  " VCR  COt  eR  (11 


“.it,  (, 

oR  RPC0  V 


cos  0„  - 2 \ esc  0 - V I 

R RT  ) R CRj 


of  . 0 / Q R \ CR 

g— = 2 esc  e cos  6 - — - — 

C0  \ T / C0 

The  chain  rule  gives  the  partials  we  need  in  terms  of  these: 


CR  3f 


30„  9v  90 

R C0  R 


3v  9v  Q 

CR  _ _9f  9f C0 

9r  9r  9v  q 9r 

Co 


Thus,  given  0V  0/9r  and  9v  Q/36  , the  required  partials  of  V can  be 
Co  Co  R CR 

determined . 

(It  may  be  of  interest  here  to  note  that  9v  /9t  , though  not  required 

CR  r 

to  compute  Q,  can  also  be  computed  easily  from 

8vcr  3f  9vc0 
9tf  = 3VC0  9tf  _ 

^C  9VCR  - 3VC0  - 

Thus,  e.g.,  -r = -r U + U can  be  found  via  these  equations.) 

dt  ^ dt^  1 dt  ^ ^ 
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We  now  turn  to  the  problem  of  finding  Bv^q/Sr  and  9Vcq/90r.  Recall 


that  av^/at^  can  be  computed  directly  once  Vcq  is  determined. 


Moreover , 


if  two  independent  equations  relating  9vcq/9r,  9vcq/90r,  and  Can 


be  found,  it  may  be  possible  to  solve  for  9V  q/9r  and  9v  Q/90„  in  terms  of 

Co  Ct)  R 


Sv^g/St^..  One  such  equation  can  be  obtained  as  follows. 


Recall  that  V 


Uj  and  VQ  = V • U2, 


while  t is  the  current  time 


(not  time-to-go) . Assume  that  V is  in  the  trajectory  plane  (so  that  Vg  is 


the  horizontal  component  of  velocity)  and  the  only  force  acting  is  that  of 

->  -+ 

Then 


gravity  (but  not  that  V = v ) . 


9v 


C6 


9v 


9t 


^ 9v  ~ 90  9v  . 9t . 

C0  9r  c6  R c6  f 

9r  9t  + 90„  9t  + 9t,  9t 


(17) 


But  clearly 


9r 

9t 


50, 


v 


R'  9t 


;,tf 

9t~ 


= -1 


(18) 


Thus 


9v 


C6 


9v 


9t 


_ce„  . !!ce^e 

9r  r 90_  R 


9v 


C0 


9t . 


(19) 


R " ' f 

Next,  from  conservation  of  the  angular  momentum  Z,  we  have 


V„  R = Z is  constant. 

O 


CO) 


Differentiating  gives 

9v 


KT  R + If  V8  " 0 


(21) 


Now  if  we  set  Vc  = V,  obviously  V^q 


V0  and  VCR 


V . Also,  since  R 


will  now  follow  a path  reaching  the  target  at  the  desired  time,  V will  con- 

/C0/ 


tinue  to  equal  V^;  thus  9vfl/9t  = 9v  g’/9t.  Eliminating  9vofi/9t  between 


■c*  e 

Equations  19  and  21  then  yields 


9v„ 


!«/!! £0.v  ) + Z£» 

CR / 9t, 


ue_ 


(22) 


J 
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Since  V can  be  zero  (i.e.,  at  apogee),  it  is  not  legitimate  to  divide 
CR 

through  by  V to  obtain  a formula  for  3v  a/3r.  This  will  cause  no  problem. 

CR  Ct/ 

Equation  22  is  the  first  of  the  two  independent  equations  sought  relating 
the  partials  of  V^g.  Since  the  Q-matrix  is  symmetric,  it  would  seem  that  the 
second  relation  needed  can  easily  be  obtained  by  setting  Qx  = Q2  i.e. 

9v, 


'C0 


3r 


i(3vcr 

r\30_ 


+ v 


C0 


(23) 


Combining  this  with  Equation  14  yields  a simple  equation  relating  Sv^q/Sr 
and  3VC0/38R: 


3VC0  l/3f  3f  3VC0  . .. 

3R  ='RK  3vc6  39R  C6, 


(24) 


Unfortunately,  when  Equations  22  and  24  are  considered  as  a system  of 

linear  equations  in  Sv^g/SR  and  3vcg/30R,  the  determinant  of  the  system  is 

zero  for  the  minimum-energy  case,  which  is  certainly  a case  of  interest.  (See 

Appendix  C.)  So  Equation  24  cannot  always  be  used  to  determine  3vcg/30R  and 

3Vcq/3r.  (However,  note  that  Equation  24,  unlike  Equation  22,  can  always  be 

used  to  find  3v  q/3r.)  Thus  another  such  relation  must  be  found. 

Ct) 

This  new  relation  will  be  obtained  indirectly  by  finding  two  independent 
equations  relating  3V(^0/3Rit,,  3vcq/3r,  3vnfi/30D,  and  SV^p/St^,  and  eliminating 


C0 

3V  a/3r  between  them. 

C0  T 


C0 


ce' 


The  first  of  these  can  be  obtained  by  considering  the  "reverse"  trajectory 

from  R to  R in  time  t . This  trajectory  passes  through  the  same  points  as 
* ~+ 
the  original  one,  but  with  the  time  sense  reversed.  That  is,  if  R(t) , V (t) 

s 

are  the  position  and  velocity  at  time  t on  the  original  trajectory,  and  R*(t), 
Vc*(t)  those  for  the  reverse  one,  then  R*(t)  = R(t^  - t)  for  0 _<  t _<  t^,  and 
thus  Vc*(t)  = -Vc(tf  - t) . 
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Now  define  V*R  to  be  the  radial  component  of  V*  at  Rt»  and  V*g  the 

horizontal  component  in  the  direction  of  R.  (Thus  V*n  and  V*  are  both 

CU  CR 

positive.)  Also,  let  VCRT»  vcqt  be  the  tadial  and  horizontal  components 
of  V(tf).  Clearly  VC0T  = V*Q  and  VCRT  = -V*R.  Since  VcQ  R = £ is  constant, 

VCQ  = VChT  = 6VC0  ' Where  6 = R/RT  * (25: 

Since  this  relationship  continues  to  hold  when  t and  6 are  allowed  to  vary, 

t R 

3u?e  . 3vce  3vce  . 3VC0 

3h  ‘ 3h  ’ % 36r  • 

The  relationship  between  Sv^q/Sr^  an<3  3v*0/3rt  as  not  ^uite  as  simple  since 
<5  depends  on  R^,  but  a simple  application  of  the  chain  rule  to  Equation  25 
yields 

■ = 6-1  + — v (27) 

3Rm  3r„  Rm  ce  1 1 

T T T 

Also,  the  same  argument  used  to  derive  Equation  22  can  be  applied  tc 
the  reverse  trajectory  to  show  that 


3v* 

V*  — 

CR  3Rt 


!Zie.v. 

HT  U„  vcr 


Substituting  Equations  25,  26,  and  27  into  Equation  28  yields 


3vce  1_/,!c0_!!c9  RT  3v0 

3RT  " M VCR*  39R  + VCR*  8tf 


(For  any  practical  case  V < 0,  so  V*  = -V  is  never  zero.) 

CRT  CR  CRT 

The  second  relation  involving  3vc0/3rt  can  derived  from  dimensional 
considerations.  Recall  that  VC0  = g(t^,  R,  Rt,  0r,  K2).  The  value  of  g does 
not  depend  on  any  physical  constants  other  than  the  arguments  themselves, 
so  that  this  function  can  be  regarded  as  a purely  numerical  relation  among 


r 


the  variables  vcq»  t^,  0 , R,  RT,  K2 . Thus,  if  the  numerical  values  of 

tf,  R,  R^,  and  K2  change  because  of  a change  in  the  physical  units  (i.e., 

units  of  distance  and  time) , the  change  in  the  value  of  will  be  the 

Co 

same  as  if  the  physical  conditions  had  changed  by  the  same  amount.  E.g., 
changing  the  time  unit  must  have  the  same  effect  on  the  numerical  value 
of  Vcg  as  the  corresponding  variations  in  the  actual  time  of  flight  (tf) 
and  the  strength  of  the  gravity  field  (K2). 

Now  suppose  we  choose  fixed  values  t , 0 , R , R , K2 , and  V „ = 

fo  R0  o To  o Co0 

g(t,_.  R /•••)  and  change  the  distance  unit  by  a factor  1/a,  so  that  all 
to  o 

distances  are  multiplied  by  a factor  a.  Since  K2  has  dimensions  of 
(distance) 3/ (time) 2 , its  value  will  be  multiplied  by  a3 . g has  dimensions 
of  distance/time,  so  its  value  will  be  multiplied  by  a.  Thus 

avce„  ' “V  “V  V C‘’K«)  . 1301 

If  we  now  change  the  time  unit  by  a factor  1/B,  similar  reasoning  yields 


?VC0o  = g(3tfo'  °V  afV  V 

The  value  of  K2  must  remain  constant,  since  otherwise  9Vcq/9k2  will 


(31) 


appear  in  the  final  result.  This  can  be  arranged  by  setting  6 = ct3/ 2 , giving 


a-1/2  v 


C0( 


- <.(a’/z 


t , aR  , aR  , 0_  , K2  \ 
to  o To  Ro  o / • 


(32) 


Differentiating  with  respect  to  a gives 


- \ a 3^2  v „ = ~ g(a3^2  t , aR  , aR  , 6 , k!) 

2 C0o  da  fo  o To  Ro  o / 


= |al/2  t 


fo  9t, 


+ R 


4*  + R. 


0 9r 


T,  3^  _ 


(33) 


where  the  derivatives  are  evaluated  at  the  point  (a3/2  t,  , aR  , aR  , 0 , K2) 

fo  o To  Ro  o 


10 


An  expression  relating  the  derivatives  at  the  point  (t,  , R. , R n , 0 


T0'  R° 1 


Kg ) can  be  obtained  by  setting  a = 1: 


_ I v = 1 t ^3_+R^.+  R 

2 c6o  2 f0  3tf  o 3r  To  3rt 


(34) 


Since  the  derivatives  here  are  evaluated  at  the  original  point,  the  zero 
subscripts  may  be  dropped  and  the  convention  of  denoting  3g/9x  by  3Vcq/3x 
(x  = t^,  R,  etc.)  may  be  resumed.  Thus  Equation  34  can  be  rewritten 


3v 


C0 


3v 


C0 


3r 


— (i  v + — t 

Rt  1 2 C0  2 3tf 


3v 


+ R 


C0 


3r 


(35) 


Eliminating  ^Vcg/^RT  between  Equations  29  and  35  yields 


30.  2 C0  |v__*  ■ 2 "flat 


t i t \#,C« 


CR 


CR 


_ 3vce  „ 
* R — ’ 0 


(36) 


3v e 

Multiplying  through  by  V and  replacing  V -3  ■ ■ by  the  riqht  side  of 

CR  CR  oR 

Equation  22  gives 


V \ 3 V _ 

1 . CR  I C6  1 

^ei1  + 5 VCR*  ) 30„  “ 2 VC0  VCR 


R + V 


T lit. 


cr1vcr*  2 


3'V 


CG 


3t, 


(37) 


It  is  shown  in  Appendix  C that  the  coefficient  of  3v  /30r  in  this 

equation  is  never  zero  for  cases  of  interest.  This  equation  can  therefore 

be  used  to  find  3v  n/30^,  in  terms  of  3v  n/3t,.  Of  course,  V*  must  first 
C0  R C0  f CR 

be  computed.  The  simplest  way  to  do  this  seems  to  be  to  apply  Equation  1 
to  the  reverse  trajectory.  This  gives 


V*  = 
CR 


cos  eR  - — lv*e  + 


1 - cos  0. 


"T  C0  \ 

Applying  Equation  25  gives  the  more  convenient  form 

,2 


V* 

CR 


(6  C°S  6R  ‘ 1)VCQ  + “ C°S  9r) 


/sin  6, 


/sin  9, 


(38) 


(39) 
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The  derivation  is  now  complete.  The  suggested  algorithm  for  computing 

Q may  be  summarized  as  follows.  Vcq  and  should  be  available  from 

the  V computation.  V and  V*  are  found  using  Equations  1 and  39  respec- 
C CR  CR 

tively.  Next  9Vcg/98R  is  computed  from  Equation  37.  The  partials  of  f are 

found  from  Equations  11,  12,  and  13,  and  9v  /90  from  Equation  14.  Now 

CR  R 

9v  „/3r  can  be  computed  from  Equation  24,  and  3V  _/3r  from  Equation  15. 

CU  CR 

Equation  10  then  gives  Q in  the  local  coordinate  frame.  Finally,  the  sym- 
metry of  Q can  be  exploited  to  save  time  in  transforming  it  to  the  desired 
frame . 


COMPUTATION  OF  Vc , VCQ,  VCR,  and  ^cQ/^{ 

First,  an  algorithm  will  be  derived  for  computing  the  value  of  t 

corresponding  to  a given  value  of  V _ for  specified  initial  and  target 

Co 

vectors.  This  derivation  essentially  follows  that  given  originally  by 
T.  Alexander.  Some  additional  terminology  must  be  introduced  here: 

e = eccentricity  of  this  trajectory  ellipse 

l = specific  angular  momentum 

0 = true  anomaly 

In  this  section,  motion  is  assumed  to  be  along  the  trajectory  ellipse. 

As  usual,  V„  and  V„  denote  the  horizontal  and  vertical  components  of  velocity, 
u R 

respectively.  However,  the  symbols  and  V will  be  reserved  for  the 

Ct)  CR 

initial  values  of  these  variables.  Initial  values  of  other  quantities  will 
be  denoted  by  zero  subscripts. 

Finally,  we  define  0„  = 0 - 9 . 0 will  be  taken  to  be  positive  in  the 

Do 

direction  of  motion,  so  that  0D  > 0 for  t > 0. 

The  usual  formula  for  R is: 

R - ■ -Q-  (Al) 

1 + e cos0 

Since  i = RVfi  = R2$,  differentiating  this  equation  gives  immediately: 


VR  = ~ 6 Sln8 


Next,  Equation  (Al)  can  be  rewritten  in  the  form: 

K2R[1  + (e  cos©  )cos0  - (e  sin0o)sin0  ] = S.2  = R„V2 


Solving  Equation  A1  for  e cos0o  and  Equation  A2  for  e sin0o  and  substi- 
tuting the  results  into  Equation  A3  yields: 


R[K2  + (R0V2fi  - K2)cos0n  - ROVc0VCRsin0n]  = R2  V2ft  - 


C0 


(A4) 


Now  at  the  target  R = R1J,  and  0Q  = 0R.  Substituting  these  values  into 
Equation  A4  and  solving  for  roVcqVcR  gives: 


R.V  „V 
0 C0  CR 


r ^)R.vJe*K2|1-c°sVl''5i"V  <A5> 

Substituting  the  expression  on  the  right  for  roVc0Vcr  2-11  E<luat^on  A4 


|cos0 


and  setting  0 = R0v£g»  we  have: 


Rl^R^sinO^  + RTsin0R(6  - K2)cos0D  - {r^B  - K2)cos0R 


- Ro0  + RTX2}sin0Dl  = 0Ro  RTsin0 


R * 


(A6) 


This  equation,  unlike  Equation  A4,  does  not  involve  V , and  thus 
gives  R explicitly  as  a function  of  0D  and  the  boundary  conditions  R0,  R^, , 

and  VC0* 

Now  l = R0vcq  = R2(d0D/dt),  or  RoVc0dt  = R2d0D,  so  that: 


/0R 

RZ(9DlljeD  . 


(A7) 


It  remains  only  to  evaluate  the  integral  in  this  equation  and  simplify 
the  resulting  expression  for  t^  as  much  as  possible.  For  this  purpose, 
define: 

A = K2RTsin0R 

B = R_(B  - K2 ) sin0 
T R 


A-2 


1 


i. 


#7 


C = 0<RO  - R.J,  coseR)  - k2r^,  (1  - cos0R) 


D = (RQ0)3/2 (RTsin0R): 


It  can  be  verified  easily  that: 


/R 

d0 

(A  + B cos©^  + 


C sin@D)  ■ 


-B  sin0„  + C cosO,. 

R R c 

F (A  + B cos0  + C sin0  ) F(A  + B) 

K R 


2A  . _i T F1/2  H 
F3/2  tan  F + C (C  + 


tan^0R 

H tan^0R) 


where 


F = A2  - B2  - C2, 


H = A - B. 


An  efficient  algorithm  for  evaluating  this  expression  for  t^  is  given 
below.  Although  a detailed  derivation  will  not  be  given,  the  following 
key  relations  were  used: 

A + B cos0„  + C sin0„  = R.  0 sin0„ 

R R o R 

H tanJj0R  = Rt(2Kz  - 3)  (1  - cos0R)  (A9) 

F + C(C  + H tan^0R)  = Rt(2K2  - 0)Rt  0 sin20R  + C(1  - cos0R)  . 

It  is  sometimes  useful  to  know  for  what  values  of  V_Q  a solution  exists. 

cy 

The  acceptable  values  are  those  on  the  open  interval  (vcgm^n*  Vc0max^ ' E<Jua" 

tions  for  V Q . and  V Q are  also  given  below. 

Cvinin  CWiricix 

The  algorithm  is  as  follows.  The  value  of  Vcg  used  is  denoted  V*g  and 
the  corresponding  value  of  t^  by  t*.  Also,  R0  is  denoted  simply  by  R,  and 
Cj  = cos0  , C 2 = sin0  , and  u = 1/K2 . 


tefc.  . . I TMiTir 


Compute : 


SI  = 1/C2 
bl  = R/RT 


d2  = 1 + bjtbj  - 2c j ) 


bs  = 1 + c, 


dl  = ^Ib5R 


bg  = i + b 

iVC0min'  VC0max  = C2/ Tdi  (b9  1 


b = b b 

6 15 


- d 

9 2 


M 


b = yTb  Rf, 

8 I 6 T 

(*)  b2  = Slv*0 

d3  = dlb2 
b7  " d3b2 


fa10  = b9  " b7d2 


bll  = V., 


b12  bll  + b7b10 


If  b,  < 0,  V*  t (V_Q_ . _ r V, 


) , so  there  is  no  solution 


12  _ *ce  " ' COrnin'  C0max  ' 

for  this  value  of  Vc0.  Otherwise: 

b.,  = Sb~ 

13  12 

d4  = bll/b13 

d5  = tt/2  - arctan  (d ^ ) (This  is  the  standard  two-quadrant  arctangent.) 
bl 4 = b6b7d5 

d6  = b 1 2b 1 3 

= bj  0bi 3 + 2b14 


= b2b8d7/d6 


V-  9 


This  completes  the  algorithm  for  evaluating  tf  from  Equation  A8.  This 
equation  can  be  differentiated  with  respect  to  V^q  to  obtain  an  expression 
for  3tf/8vcQ,  and  thus  of  = (3tf/3Vc0) -1 . Again  without  going 

through  the  algebraic  details,  an  algorithm  for  evaluating  this  expression 
is  given  below.  The  resulting  value  of  is  denoted  Qv^g/St^).* 

Compute : 


c 

= 2d,  S 

J | 

7 

3 I 

' 1 

C 1 2 

= 2C7bi  0 

C 1 3 

= c12/(2b13) 

'' 

C 

= b)cd  +bd(dc  -cb  ) / [b  (1  + d2 ) f 

14 
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! 

‘8W3tf>* 

d6/|be[Sxd7  + b2^b10C13  + 2C14  ~ b13C7d2‘l] 

! 

- t*(c 

1 2bZ  3 + bl  2C1  3 

|' 

• 

-r-  ->■  o 

In  general,  the  values  given  initially  are  R,  R^,  K , and  tf,  while  the 

values  required  are  V and  possibly  V_Q,  V , and  3v  /3t  . The  algorithms 

C Cu  CR  Cu  t 

above  give  tf  and  3Vc0/3tf  explicitly  in  terms  of  R, 

VCQ.  R,  Rt,  cos0R,  and  sin0R  can  be  found  immediately  from  the  given  values, 

while  V can  be  determined  once  V„_  and  V are  known.  Also,  V _ can  be  com-  I 

C C0  CR  CR 

puted  for  any  given  value  of  from  Equation  1 which  can  be  written  in  the 

form:  I 

VCR  = V(b8b2}  + b2(cl  ~ V . 

NOTE:  Since  the  value  of  V„Q  finally  accepted  may  not  be  the  last  value  of 

Co 

V*Q,  = S_V_Q  must  in  general  be  recomputed  before  V is  evaluated. 

Cu  2 I Cu  CR 


Rt,  COS0R,  sin0R,  and 


A-5 


“ * ' 


Thus  the  problem  is  reduced  to  that  of  finding  the  value  of  V^g  cor- 
responding to  the  given  value  of  t^.  The  required  value  must  be  found  by 
iteration.  Since  t^  is  a well-behaved  function  of  V^g  and  3Vcg/3t^  is 
available,  the  Newton- Raphson  method  is  an  obvious  possibility.  However, 
since  A = (Sv^g/St^ ) t^/sin6R  is  much  more  nearly  constant  than  3Vcg/3tf, 
more  rapid  convergence  can  be  obtained  by  choosing 


vJe  * <sW8tf)*lt£/tfHtf 


tp  , 


where 


t^  is  the  given  value  of  time  of  flight, 
as  the  next  estimate  of  Vcg  on  each  iteration. 

Only  the  formulas  from  (*)  on  need  to  be  reevaluated  on  each  iteration. 

If  no  better  initial  estimate  of  V^g  is  available,  the  value 

vce  ■ vc6„in  + *r  sln6R/t£ 
is  usually  satisfactory. 

Lastly,  if  the  final  value  of  Vcg  is  not  the  last  value  of  V*g  used  in 
the  algorithm,  a better  estimate  of  the  corresponding  value  of  ^vcg/3tf  than 
O^g/atf)*  is 

3vc6/3tf  - Ovc(,/3tf)*(t*/tf)J  . 


DERIVATION  OF  EQUATION  FOR  V 


The  importance  of  the  Q-matrix  stems  in  large  part  from  its  occurrence 

4- 

in  the  well-known  equation  for  V : 

g 


-i  -+ 

V = -A 

g t 


Qv 


(Bl) 


where  At  is  the  sensed  (or  gravity-free)  acceleration;  i.e.,  the  acceleration 
produced  by  all  forces  except  gravity.  This  result  is  established  by  the 
following  simple  argument. 

First,  with  the  conventions  explained  in  the  Introduction,  clearly 


4 - 3VC 

vc  = 2V  ‘ 9T 


(B2) 


On  the  other  hand,  V = A^  + G by  definition  of  At.  An  equation  relating 
G to  V and  its  derivatives  can  be  obtained  most  simply  by  considering  the 

-V  ->■ 

"target"  trajectory;  i.e.,  the  free-fall  trajectory  for  which  V = V initially. 

As  pointed  out  in  the  Introduction,  on  this  trajectory  V = V for  all  t,  so 
-V  4-  ->■  4-  L 

that  = V = G.  Combining  this  with  the  expression  for  given  by  Equa- 
tion B2  yields 


9v 

-y  -y  p 

G = QVc  - 9^7 


(B3) 


Now  recall  that  = V^tt^,  R,  R^,  K2),  so  that  neither  nor  any  of 
its  derivatives  depend  on  the  values  of  V or  A . Thus,  while  Equation  B3 

■+  ~y 

is  most  easily  derived  by  assuming  that  V = V and  A = 0,  actually  both 

~y  ~y  ^ 

sides  are  independent  of  V and  A^,  so  that  this  equation  is  valid  for  any 
values  of  these  variables.  We  can  therefore  write 


V = A + G = A +QV 
T T C 


9t„ 


(B4) 


Subtracting  Equation  B4  from  Equation  B2  yields  Equation  Bl. 


B-l 


APPENDIX  C 


NECESSITY  AND  SUFFICIENCY  OF  EQUATIONS 
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NECESSITY  AND  SUFFICIENCY  OF  EQUATIONS 


This  section  deals  with  the  suitability  of  the  algorithm  developed  in 
the  main  section  and  discusses  alternative  methods  of  computing  the  partials 
of  V 


C0" 


The  most  important  question  to  be  considered  is  whether  the  coefficient 
of  ^Vcg/^®R  i-n  Equation  37  can  ever  be  zero  in  any  case  of  interest.1  Since 


Vcg  is  always  positive,  this  is  equivalent  to  asking  whether  it  is  possible 


that 


i + #-^-  o 


R V * 
T CR 


(Cl) 


for  such  a case.  Since  V * > 0 for  all  such  trajectories,  clearly  V must 

CR  CR 


be  negative  for  Equation  Cl  to  hold;  i.e.,  the  body  must  be  past  apogee. 


Since  V = -V  *,  Equation  Cl  is  equivalent  to 
CRT  CR 


rvcr  rtvcrt 


(C2) 


But  by  Equation  20,  in  any  case 


RVce  = Vcqt 


(C3) 


Squaring  Equations  C2  and  C3  and  adding  gives 


R?V  = R^  vi_ 

C T CT 


(C4) 


where  VCT  = |Vc(tf) |. 


By  the  well-known  vis-viva  equation. 


V‘  = 1C  (2/R  - 1/a)  , 


(C5) 


where  a is  the  semimajor  axis  of  the  conic.  Since  a is  constant  and  Equa- 


tion C5  must  hold  at  both  t = 0 and  t = t^,  this  gives  expressions  for 


V and  VCT  which  can  be  substituted  into  Equation  C4  to  yield 


(R  - Rt) [2a  - (R  + RT)]  = 0 . 


(C6) 


1 The  term  "case  of  interest"  refers  here  to  conditions  which  might  arise 
during  powered  flight  for  a ballistic  missile. 


C-l 


! 

A 


The  only  time  R can  equal  past  apogee  is  at  impact.  At  this  point 
calculations  are  no  longer  required,  so  this  case  can  be  ignored.  Otherwise 

the  expression  in  brackets  must  be  zero.  Thus  we  must  have  R + R = 2a,  or 

1 T 
a = — (R  + R^) . We  will  show  that  computation  of  3Vcq/30r  will  never  be  re- 
quired past  apogee  for  such  a trajectory. 


The  sum  of  the  distances  from  any  point  on  an  ellipse  to  the  two  foci 

is  2a.  Thus  the  distance  from  R to  the  empty  focus  must  be  R . This  is 

illustrated  in  Figure  1.  But  between  apogee  and  impact  R > R^,  so  that  RT 

must  lie  on  the  same  side  of  the  minor  axis  as  0,  the  center  of  the  earth. 

- 

Now  let  Rq  be  the  other  point  on  the  ellipse  whose  distance  from  0 is  R . 

The  time  required  to  go  from  R0  to  R^  is,  by  Kepler's  second  law,  (Ag/AE)P 

where  A is  the  shaded  area  in  the  diagram,  A is  the  total  area  of  the 

ellipse,  and  P is  the  period.  Clearly  A /A  > 1/2,  and  a well-known  formula 

S E 

for  P is 


P = 2JL  a3/2 


(C7) 


Now  since  R > R^,  a = — (R  + R^)  > R^,  and  since  R^  must  be  essentially 

on  the  earth's  surface  we  have  a > 2 x 107  ft.  Also,  K2  < 1.5  x 1016  ft3/sec2, 

-y  ~y 

so  P > 4000  sec.  Thus  the  time  required  to  go  from  R0  to  R^  is  at  least 
2000  sec,  and  the  time  to  reach  apogee  is  half  that,  or  at  least  1000  sec. 

Since  the  actual  launch  point  is  also  on  the  earth,  the  time  required  to 
reach  apogee  on  the  actual  trajectory  is  clearly  greater  than  that  required 
to  reach  apogee  from  RQ  on  the  ellipse.  But  this  formula  will  never  be 
needed  as  late  as  1000  seconds  after  launch.  Since  Equation  Cl  can  only 
be  satisfied  past  apogee,  this  shows  that  it  cannot  be  satisfied  in  any 
case  of  interest. 


This  shows  that  Equation  37  can  be  used  to  compute  9vcg/99R  in  a 9eneral 
algorithm.  However,  the  question  of  whether  Equation  24  can  be  used  to  pro- 
duce a (possibly  simpler)  equation  for  this  purpose  should  also  be  considered. 


C-2 


TRAJECTORY  ELLIPSE 


0 = CENTER  OF  EARTH 
E = EMPTY  FOCUS  OF  TRAJECTORY  ELLIPSE 


OR+RE  « 2a  = R+Rt 
OR  = R 
•••  RE  « Ry 

RE  < OR 


Figure  1.  Case  Which  Satisfies  Equation  Cl 


If  Bv^q/Sr  is  eliminated  between  Equations  22  and  24,  the  result  is 


V + V JLL. 

ce  CR  3vce 


v !!_  + R 

cr  aeR  K at 


Thus  the  question  is  whether  the  coefficient  of  ^vcq/^®r  in  this 
equation  can  ever  be  zero;  i.e.,  whether  the  equation 

vce  + v«  #7  ’ 0 <«> 

Cu 

can  be  satisfied  in  any  case  of  interest.  Since  V = V U + V « U-? , clearly 

C CR  1 C0  2 1 

Vc  = VC6  + VCR  . (cl0) 

As  explained  in  the  DERIVATION  section,  Equation  2 defines  the  function 
f<Vc6'  R'  V ®r*  Thus,  by  Equation  CIO,  can  also  be  regarded  as 

a function  of  these  parameters.  Differentiating  with  respect  to  yields 

8vc  3f 

vc^=vce  + vcR3^  . (cll) 

Thus,  Equation  C9  is  satisfied  if  and  only  if  = 0.  This 

will  happen  for  any  value  of  V Q at  which  V_  attains  a minimum  (or  maximum) 

Cb  C 

value;  in  particular,  for  the  value  of  V^g  corresponding  to  the  minimum 
energy  trajectory.  This  is  certainly  a case  of  interest;  in  fact,  a tra- 
jectory close  to  this  one  is  often  chosen  deliberately  to  minimize  the 
energy  needed  to  reach  the  target. 

(Clearly  the  minimum  energy  trajectory  must  be  an  ellipse.  Each  value 

of  V . on  the  open  interval  (V  a , V n ) corresponds  to  a unique  ellip- 
Co  ^ CUmax 

tical  path  from  R to  R , and  vice-versa,  so  the  fact  that  V„  is  a minimum 

T C 

does  in  fact  show  that  = 0.) 

We  have  shown  that  Equation  C8  in  itself  is  not  suitable  for  use  in  a 
general-purpose  algorithm  for  computing  Q.  Still,  it  might  be  possible  to 
use  Equation  C8  together  with  Equation  37  to  determine  Dv^g/SOp  and  ^Vc0/,^tf 
simultaneously.  This  would  certainly  be  useful  in  some  circumstances.  For 


example,  in  some  applications  a good  estimate  of  v^g  may  be  available  without 
using  the  usual  Vc  iteration;  in  this  case  3Vcg/3tf  is  not  automatically 
available. 

Equations  28  and  C8  can  be  combined  into  the  matrix  equation 


''cel1  + s V". 

Cft 


3f 

V + V -x — ~ 

C0  CR  3V 


R + V 


CRl  V 


ce 


r 3vc6_ 

30 

R 

3vce 

. 

L9tf  . 

— V V 
2 C0  CR 


-V 


3f 
CR  96 


R J 


(C12) 


This  system  can  be  solved  if  the  determinant  of  the  matrix  on  the  left 
is  not  zero.  Thus  the  question  to  be  considered  is  whether  the  equation 


RVcel1  1 s v 


CR 


CR 


= I V + V 

1 ce  CR  3v 


9f 


C0 


R + V 


CR\V 


(C13) 


can  be  satisfied  for  any  cases  of  interest.  This  equation  can  be  satisfied, 
but  it  is  not  known  whether  any  situation  in  which  it  holds  could  actually 
occur  in  practice. 

A further  question  is  whether  still  another  relationship  among  Sv^q/Sr, 
9Vcq/90r,  and  3Vcg/3tf  independent  of  those  derived  here  might  not  be  found. 
Several  fairly  obvious  ways  to  derive  such  relations  besides  those  described 
in  this  paper  can  be  applied,  but  they  all  lead  to  relations  which  follow 
from  those  already  obtained.  On  reflection,  it  can  be  seen  that  this  was  to 
be  expected.  Methods  along  the  lines  of  those  given  here  will  lead  to  linear 
relations  among  these  derivatives.  Another  such  relation  would  make  it  possible 
to  eliminate  all  the  derivatives  of  V^g  to  give  an  algebraic  expression  re- 


lating tf,  Vcg,  VCR,  R,  Rt,  K , sin  8r,  and  cos  0R. 


The  occurrences  of  V 


CR 


could  then  be  eliminated  by  applying  Equation  1.  The  result  would  be  an 


algebraic  relation  between  t^  and  V^g  in  terms  of  the  other  parameters  listed. 


f 

[ 

! But  this  impossible,  because  it  would  lead  to  the  result,  for  example,  that 

[ and  arctan  d4  are  also  related  in  this  way.  (See  Appendix  A.) 

Nevertheless,  there  are  some  problems  to  be  resolved.  One  is  to  find  a 

t 

useful  characterization  of  the  conditions  under  which  Equation  C13  is  satisfied. 
Another  is  to  determine  whether  a relation  between  3Vcq/90r  and  3VcQ/3tf  can 
be  found  which  does  not  involve  tf  explicitly,  and  in  which  the  coefficient 
of  3Vc0/30r  is  never  zero.  (This  would  be  useful  in  certain  applicatons) . 

, Similarly,  alternative  expressions  relating  the  derivatives  of  VqQ , though 

equivalent  to  those  used  here,  might  be  of  such  a nature  as  to  permit  simul- 
taneous solution  for  all  cases  of  interest  (or  ideally,  for  all  cases) . These 
questions  are  appropriate  subjects  for  further  research. 
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